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The response of complex networks to perturbations is of utmost importance 
in cireas as diverse as ecosystem management, emergency response, and cell 
reprogramming. A fundamental property of networks is that the perturbation 
of one node can affect other nodes, in a process that may cause the entire or 
substantial part of the system to change behavior and possibly collapse. Recent 
reseau"ch in metabolic and food-web networks has demonstrated the concept 
that network damage caused by external perturbations can often be mitigated 
or reversed by the application of compensatory perturbations. Compensatory 
perturbations are constrained to be physically admissible and amenable to 
implementation on the network. However, the systematic identification of 
compensatory perturbations that conform to these constraints remains an open 
problem. Here, we present a method to construct compensatory perturbations 
that can control the fate of general networks under such constraints. Our 
approach accounts for the full nonlinecir behavior of real complex networks 
and can bring the system to a desirable target state even when this state is 
not directly accessible. Applications to genetic networks show that compen- 
satory perturbations are effective even when limited to a small fraction of all 
nodes in the network and that they are far more effective when limited to the 
highest- degree nodes. The approach is conceptually simple and computationally 
efficient, making it suitable for the rescue, control, and reprogramming of Icirge 
complex networks in veo-ious domains. 



Introduction 



Complex systems such as power grids, cellular networks, and food webs are often modeled 
as networks of dynamical units. In such complex networks, a certain incidence of pertur- 
bations and the consequent impairment of the function of individual units — whether genes, 
power stations, or species — are largely unavoidable in realistic situations. While local per- 
turbations may rarely disrupt a complex system, they can propagate through the network 
as the system accommodates to a new equilibrium. This in turn often leads to system-wide 
re-configurations that can manifest themselves as genetic diseases [ll[2], power outages |3lll], 
extinction cascades [HI E], traffic congestions [U E], and other forms of large-scale failures 

pun]. 

A fundamental characteristic of large complex networks, both natural and man-made, is 
that they operate in a decentralized way. On the other hand, such networks have either 
evolved or been engineered to inhabit stable states in which they perform their functions 
efficiently. The existence of stable states indicates that arbitrary initial conditions converge 
to a relatively small number of persistent states. Such states are generally not unique and can 
change in the presence of large perturbations. Because complex networks are decentralized, 
upon perturbation the system can spontaneously go to a state that is less efficient than 
others available. For example, a damaged power grid undergoing a large blackout may still 
have other stable states in which no blackout would occur, but the perturbed system does 
not converge to those states. We suggest that many large-scale failures are determined by 
the convergence of the network to a "bad" state rather than by the unavailability of "good" 
states. 

Here we explore the hypothesis that one can design physically admissible compensatory 
perturbations to direct a network to a desirable state even when it would spontaneously go 
to an undesirable ( "bad" ) state. An important precedent comes from the study of food- web 
networks, where perturbations caused by invasive species or climate change can lead to the 
subsequent extinction of numerous species. Recent research predicts that a significant frac- 
tion of these extinctions can be prevented by the targeted suppression of specific species in 
the system [TT]. Another precedent comes from the study of metabolic networks of single-cell 
organisms, where perturbations caused by genetic or epigenetic defects can lead to nonviable 
strains, which are unable to reproduce. The knockdown or knockout of specific genes has 
been predicted to mitigate the consequence of such defects and often recover the ability of 
the strains to grow [12]. These findings have analogues in power grids, where perturba- 
tions caused by equipment malfunction or operational errors can lead to large blackouts, 
but the appropriate shedding of power can substantially reduce subsequent failures [I3|, . 
Therefore, the concept underlying our hypothesis is supported by recent research on physical 
[mdS], biological [121 US], and ecological networks [H]. In these studies, the identification 
of compensatory perturbations relied either on simple models of the real system [HI [15] 
or on using heuristic arguments [TTl [12] . The question we pose is whether compensatory 
perturbations can be systematically identified for a general network of dynamical units. 

Our solution to this problem is based on the observation that associated with each de- 
sirable state there is a region of initial conditions that converge to it — the so-called "basin 
of attraction" of that state. Given a network that is at (or will approach) an undesirable 
state, the conceptual problem is thus reduced to the identification of a perturbation in the 
state of the system that can bring it to the attraction basin of the desired stable state (the 
target state). Once there, the system will evolve spontaneously to the target. However, this 
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perturbation must be physically admissible and is therefore subject to constraints — in the ex- 
amples above, certain genes can be down-regulated but not over-expressed, the populations 
of certain species can only be reduced, and changes in power flow are limited by capacity 
and generation capability. Under such constraints, the identification of a point within the 
basin of attraction of the target state is a highly nontrivial task. This is so because complex 
networks are high-dimensional dynamical systems and it is generally impossible to identify 
the boundaries of a basin of attraction in more than just a handful of dimensions. 

The approach introduced here overcomes this difficulty and is able to identify com- 
pensatory perturbations in the absence of any a priori information about the location of 
the attraction basin of the target state. This is done systematically, without resorting to 
trial-and-error. The approach is effective even when only exploiting resources available 
in the network, which usually forbids bringing the system directly to the desired stable 
state. This often limits the search to perturbations that are locally deleterious, such as the 
temporary inactivation of a node. Interestingly, applications of our method show that such 
locally deleterious perturbations can be tailored to lead to globally advantageous effects. 
The approach is general and can be used to both control a network in the imminence of a 
failure and to reprogram a network even in the absence of any impending failure. 

Results 

Figure [T]A.-C illustrates the problem that we intend to address. The dynamics of a network 
is best studied in the state space, where we can follow the time evolution of individual 
trajectories and characterize the stable states of the whole system. Figure represents 
a network that would spontaneously go to an undesirable state, and that we would like 
to bring to a desired stable state by perturbing at most three of its nodes (highlighted). 
Using X = (xi,...,x„) to represent the dimensions of the state space. Fig. [l|3 shows how 
this perturbation applied at a certain time to 5 changing the state of the system from xq to 
Xq, would lead to an orbit that asymptotically goes to the target state. As an additional 
constraint, assume that the activity of the nodes can only be reduced (not increased) by this 
perturbation. Then, in the subspace corresponding to the nodes that can be perturbed, the 
set of points X that can be reached by eligible compensatory perturbations forms a cubic 
region, as shown in Fig. [ip. The target state itself is outside this region, meaning that it 
cannot be directly reached by any eligible perturbation. However, the basin of attraction 
of the target state may have points inside the region of eligible perturbations (Fig. [lf3), 
in which case the target state can be reached through perturbations that bring the system 
to one of these points; once there, the system will spontaneously evolve towards the target 
state. This scenario is representative of the conditions under which it is important to 
identify compensatory perturbations, and leads to a very clear conclusion: a compensatory 
perturbation exists if and only if the region formed by eligible compensatory perturbations 
overlaps with the basin of attraction of the target. 

However, there is no general method to identify basins of attraction (or this possible 
overlap) in high-dimensional state spaces typical of complex networks. Indeed, despite 
significant advances, numerical simulations are computationally prohibitive and analytical 
methods, such as those based on Lyapunov stability theorems, offer only rather conservative 
estimates and are not yet sufficiently developed to be used in this context |T71 [T8] . 
Accordingly, our approach does not assume any information about the location of the 
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attraction basins and addresses a problem that, as further explained below, cannot be 
solved by existing methods from control theory. 



Control procedure for networks 

The dynamics of a complex network can often be represented by a set of coupled ordinary 
differential equations. We thus consider an A^-node network whose n-dimensional dynamical 
state X is governed by 



The example scenario we envision is the one in which this network has been perturbed at 
a time prior to to, driving it to a state xq = x(to) in the attraction basin Q{'Ku) of an 
undesirable state x^. The compensatory perturbation that we seek to identify consists of a 
judiciously-chosen perturbation xq — )■ Xq to be implemented at time to, where Xq belongs to 
the basin of attraction fi(x*) of a desirable state x*. For simplicity, we assume that x^ and 
X* are fixed points, although the approach we develop extends to other types of attractors. 
In the absence of any constraints it is always possible to perturb Xq such that Xq = x*. 
However, as illustrated above, usually only constrained compensatory perturbations are 
allowed in real networks. These constraints encode practical considerations and often take 
the form of mandating no modification to certain nodes, while limiting the direction and 
extent of the changes in others. The latter is a consequence of the relative ease of removing 
versus adding resources to the system. We thus assume that the constraints on the eligible 
perturbations can be represented by the vectorial expressions 



where both the equality and the inequality are assumed to apply to each component. In a 
food web, for example, admissible compensatory perturbations may only reduce the popula- 
tions of certain unendangered species, while leaving those of the endangered ones unchanged. 
For concreteness, we focus on constraints of the form Xqj < Xqj for accessible nodes, and 
x'qj = Xoj for the other nodes, although the approach we develop accommodates general, 
nonlinear constraints. 

We propose to construct compensatory perturbations iteratively from small perturba- 
tions, as shown in Fig. [ip-E. Given a dynamical system in the form ([T| and an initial state 
Xq at time to, a small perturbation 5xo evolves to 5x(t) = M(xo,t) ■ ^xq at time t. The 
matrix M(xo,t) is the solution of the variational equation dM./dt = DF{x.) ■ M subject to 
the initial condition M(xo,0) = 1, and is generally invertible [19]. Given the time tc of 
closest approach of the perturbed orbit to the target x*, we can in principle use the inverse 
transformation. 



to determine the perturbation ^xq to the initial condition xq that, among the admissible 
perturbations satisfying |5xo| < ei, will render x(tc) + 5x(tc) closest closest to x* (Fig. 

This is a good approximation for small ei, and hence for small changes in the initial 
conditions; the error can in fact be calculated from the matrix of second derivatives and 
verified numerically by direct integration. Large perturbations can then be built up by 
iterating the process: every time (5xo is calculated, the current initial state, Xq, is updated 
to Xq + 5xo, and a new Sxq is calculated starting from the new initial condition (Fig. fife). 




(1) 




(2) 




(3) 
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The problem of identifying a perturbation (5xo that incrementally moves the orbit toward 
the target under the given constraints is then cast as a constrained optimization problem (see 
Methods). Once found, the optimal incremental perturbation is applied to the current initial 
condition, Xq — )■ Xq + (5xo, giving the initial state for the next iteration. At this point, we test 
whether the new state lies in the target's basin of attraction by integrating the system ([T]) 
over a long time r. If the resulting orbit reaches a small ball of radius k around the target, 
we declare success and terminate the procedure. Otherwise, we identify the closest approach 
point in a time window to < t < to + T, where T is a time limit determined beforehand, and 
repeat the procedure. Now, it may be the case that no compensatory perturbation can be 
found, e.g., if the feasible region X does not intersect the target basin fi(x*). To account 
for this, we automatically terminate our search if the system is not controlled within a fixed 
number of iterations (see Methods). 

Before considering the application to networks, we first illustrate our method using 
an example in two dimensions, where the basins of attraction (and hence the possible 
compensatory perturbations) can be easily calculated and visualized. Figure [2] shows the 
state space of the system, which has two stable states: x^ on the left and x^ on the right. 
The system is defined by the potential U{xi) = exp(— 7x^)(6a;^ + cxf + dxf) and frictional 
dissipation t], where 7 = 1, 6 = —1, c = —0.1, d = 0.5, r] = 0.1, and X2 = xi. The method is 
illustrated for two different initial states under the constraint that admissible perturbations 
have to satisfy Xq < xq, i.e., one cannot increase either variable. For the initial state in the 
basin of state x^ (Fig. [2|\.), no admissible perturbation exists that can bring the system 
directly to the target x^, on the right, since that would require increasing Xi. However, 
our iterative procedure builds an admissible perturbation vector that shifts the state of the 
system to a branch of the basin ^(xb) lying on the left of that point. Then, from that 
instant on the autonomous evolution of the system will govern the trajectory's approach 
to the target x^, on the right. This example illustrates how compensatory perturbations 
that move in a direction away from the target — the only ones available under the given 
constraints — can be effective in controlling the system, and how they are identified by our 
method. The other example shown illustrates a case in which the perturbation to an initial 
state on the right crosses an intermediate basin, of x^, before it can reach the basin of the 
target, (Fig. [2|3). The linear approximation fails at the crossing point, but convergence 
is nonetheless assured by the constraints imposed on 5xo (see Methods). 



Control of genetic networks 

We now turn to compensatory perturbations in complex networks. We focus on networks 
of diffusively coupled units, a case that has received much attention in the study of spon- 
taneous synchronization [20]. We take as a base system the genetic regulatory subnetwork 
shown in Fig. [sjA. (inset), consisting of two genes wired in a circuit. The state of the system 
is determined by the expression levels of the genes, represented by the variables Xi,X2 > 0. 
The associated dynamics obeys 



dx 1 x-i ^ S ^ „ 

—r- = Ol 7; \- Oi KiXi + fi, (4) 

dX2 X^ ^ S"" , ^ . , 

where the first two terms for each gene capture the self-excitatory and mutually inhibitory 
interactions represented in Fig. [3]A., respectively, while the final two terms represent linear 
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decay and a basal activation rate of the associated gene's expression. Models of this form 
have been used extensively to describe the transition between progenitor stem cells and 
differentiated cells [2IH23]- For a wide range of parameters, this system exhibits three stable 
states! a state (xb) characterized by comparable expression of both genes (the stem cell 
state), and two states (x^ and x^r) characterized by the dominant expression of one of 
the genes (differentiated cell states). Figure [s] and subsequent results correspond to the 
symmetric choice of parameters ai 2 = 0.5, 61,2 = 1, /ci,2 = 1, /i,2 = 0.2, S = 0.5, and 
m = 4. We assume that compensatory perturbations are limited to gene down-regulations, 
i.e., Xq < xq. Admissible compensatory perturbations do not exist (and hence cannot be 
identified) from some other initial states close to the axis in the basin of ^a,c having X(7 ,4 as 
target, indicating that the expression of that gene is too low to be rescued by compensatory 
perturbations. But x^ can be reached even in these cases, and thus so can x^^ ^; if we let 
the system evolve after a perturbation into ^{'Xb) and then perturb it again into fli^c^A), 
meaning that we need the conversion from the differentiated state to the stem cell state 
before going to the other differentiated state. 

We construct large genetic networks by coupling multiple copies of the two-gene system 
described above. Such networks may represent cells in a tissue or culture coupled by means 
of factors exchanged through their microenvironment or medium. Specifically, we assume 
that each copy of this genetic system can be treated as a node of the larger network. The 
dynamics of a network consisting of such systems is then governed by 

^ = f (xi) + J- ^ Aij [xj- - Xi] , (6) 

* j 

where Xj = f(xj) is the vectorial form of the dynamics of node i as described by Eqs. (|4])- 
(|5|, the parameter e > is the overall coupling strength, and di is the degree (number of 
connections) of node i. The structure of the network itself is encoded in the adjacency matrix 
A = (Aij). We focus on randomly generated networks with both uniform and preferential 
attachment rules (Methods), which serve as models for homogeneous and heterogeneous 
degree distributions, respectively. We use x = (x,) to denote the state of the network, 
with x^, X.B, and x^ denoting the states in which all nodes are at state x^, x^, and x^;, 
respectively. The states x^i, x^, and xc are fixed points of the full network dynamics in 
the 2A^-dimensional state space and, by arguments of structural stability, we can conclude 
they are also stable and have qualitatively similar basins of attraction along the coordinate 
planes Xj if the coupling strength e is weak. While we focus on these three states, it follows 
from the same arguments that in this regime there are 3^—3 other stable states in the 
network. Under such conditions, compensatory perturbations between x.^, x^, and xc" are 
guaranteed to exist, and hence this class of networks can also serve as a benchmark to 
test the effectiveness and efficiency of our method in finding compensatory perturbations in 
systems with a large number of nodes. 

A general compensatory perturbation in this network is illustrated in Fig. [3^, where 
different intensities indicate different node states. Applied to the initial state x.a and target 
Xfi for e = 0.05, the method is found to be effective in 100% of the cases for the 10, 000 
networks tested, with A^ ranging from 10 to 100. Moreover, the computation time and 
number of iterations for these tests confirm that our method is also computationally efficient 
for large networks. Computation time grows approximately quadratically with A^ (Fig. 
[3p), as expected since each iteration requires the integration of 0{N'^) equations. The 
number of iterations grows as the square root of A^ (Fig. pb), in agreement with the \/N 
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scaling of |xa — ^b\ and of the distances between other invariant sets. This leads to the 
asymptotic scaling N^^"^ for the computation time, which is not onerous since the control 
of one network requires the identification of only one compensatory perturbation. This 
should be compared to the 0(exp(A^)) time that would be required to determine the basin of 
attraction at fixed resolution by exhaustive sampling of the state space. These properties are 
representative of more general conditions. The efficiency and effectiveness of our approach 
are expected to be relevant not only for large genetic networks, but also for large real 
networks in general, such as food webs and power grids, where the characteristic time to 
failure after the initial perturbation is large relative to the expected time scales to identify 
compensatory perturbations. 

In the preceding analysis, we allowed all nodes in the network to be perturbed, thus 
providing an upper bound for the computational time. In a realistic situation we may assume 
that only a subset of the network is made available to perturbations. While compensatory 
perturbations are not expected to exist under such constraints for very small e, they become 
abundant for larger coupling strengths. Figure |4] shows the success rate in directing the 
system from xq = to x* = for several network sizes and a varying number of nodes k 
in the control set, i.e., the set of nodes available to be perturbed. For e = 1, as considered in 
these case, the success rate increases monotonically with k, as expected, but it is frequently 
possible to identify an eligible compensatory perturbation using a small number of nodes 
in the network. For example, for networks of 50 nodes, approximately 40% of all randomly 
selected 10-node control sets will lead to successful control, for both homogeneous (Fig. |4]A.) 
and heterogeneous (Fig. |4|3) networks. Naturally, for the network to be controlled it suffices 
to find a single successful control set, indicating that the networks will often be rescued by a 
much smaller fraction of nodes. From a comparison between Fig. |4]A. and Fig. it follows 
that, for control sets formed by randomly selected nodes, heterogeneous networks are more 
likely to be controlled than homogeneous ones if the size of the control set is small, while 
the opposite is true if the control set is large. The analysis below will allow us to interpret 
this property of heterogeneous networks as a consequence of the possibility of having a very 
high-degree node in small control sets and a tendency to have a large number of low-degree 
nodes in large control sets. But what distinguishes a successful control set from a set that 
fails to control the network? 

Figure |5] shows that the probability of success is strongly determined by node degree. This 
is clear from the probability that a node appears in a successful control set as a function 
of its degree d (Fig. [5|\), which increases approximately linearly for both homogeneous and 
heterogeneous networks. This is even more evident in the dependence of the success rate on 
the average degree {d) of the control set (Fig. [5^). The success rate exhibits a relatively sharp 
transition from zero to one as (d) increases, which starts near the peak of the distribution 
of average degrees for random control sets (Fig. [5|3, background histograms). As shown in 
Fig. [5p for a homogeneous network, the success rate of control sets involving a given node 
correlates strongly with degree even across relatively small degree differences. This conffims 
both that compensatory perturbations can be identified for a large fraction of control sets 
and that successful control sets show an enrichment of high-degree nodes. 

Because successful control sets tend to have a significantly higher fraction of high-degree 
nodes, increased success rate can be achieved by biasing the selection of control-set nodes 
towards high degrees. This is demonstrated in Fig. |6j where the rate of success increases 
dramatically when the control set is formed by the highest-degree nodes in the network. 
For example, for the control sets formed by 20% of the nodes, the success rate approaches 
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100% for all network sizes considered, more than doubling when compared to randomly 
selected control sets. This increase is more pronounced for more heterogeneous networks, as 
illustrated for control sets limited to 10% of the nodes (Fig. [6]), which is expected from Fig. 
[5] and the availability of higher degree nodes in such networks. This enhanced controllability 
of degree-heterogeneous networks for targeted node selection is important in view of the 
prevalence of heterogeneous degree distributions in natural and engineered networks. We 
emphasize that these conclusions follow from our systematic identification of compensatory 
perturbations, and that they cannot be derived from the existing literature. 

Note that our approach is fundamentally different from those usually considered in 
control theory, both in terms of methods and applicability. Optimal control [2lj, for 
example, is based on identifying an admissible (time-dependent) control u(t) such that the 
modified system dyi/dt = F + u will optimize a given cost function. Control of chaos ^3] . 
which is used to convert an otherwise chaotic trajectory into a periodic one, is based on the 
continuous application of unconstrained small time-dependent perturbations to align the 
stable manifold of an unstable periodic orbit with the trajectory of the system. Approach 
to the relevant orbit can be facilitated by the method of targeting [26], which, like control 
of chaos itself, can only be applied to move within the same ergodic component. Here, in 
contrast, we bring the system to a different component, to a state that is already stable, and 
we do so using one (or few) finite-size perturbations, which are forecast-based rather than 
feedback-based. Our approach is suited to a wide range of processes in complex networks, 
which are often amenable to modeling but offer limited access to interventions, making the 
use of punctual interventions that benefit from the natural stability of the system far more 
desirable than other conceivable alternatives. 

Discussion 

The dynamics of large natural and man-made networks are usually highly nonlinear, 
making them complex not only with respect to their structure but also with respect to 
their dynamics. This is particularly manifest when the networks are brought away from 
equilibrium, as in the event of a strong perturbation. Nonlinearity has been the main obstacle 
to the control of such systems. Progress has been made in the development of algorithms for 
decentralized communication and coordination [27] , in the manipulation of two-state Boolean 
networks [28], in network queue control problems [29], and other complementary areas. Yet, 
the control of far-from-equilibrium networks with self-sustained dynamics, such as metabolic 
networks, power grids, and food webs, has remained largely unaddressed. Methods have 
been developed for the control of networks hypothetically governed by linear dynamics [30] . 
But although linear dynamics may approximate an orbit locally, it does not permit the 
existence of different stable states observed in real networks and does not account for basins 
of attraction and other global properties of the state space. These global properties are 
crucial because they underlie network failures and, as shown here, also provide a mechanism 
for the control of numerous networks. 

The possibility of directing a complex network to a predefined dynamical state offers 
an unprecedented opportunity to harness such highly structured systems. We have shown 
that this can be achieved under rather general conditions by systematically designing com- 
pensatory perturbations that, without requiring the computationally-prohibiting explicit 
identification of the basin boundaries, take advantage of the full basin of attraction of the 
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desired state, thus capitalizing on (rather than being obstructed by) the nonhnear nature of 
the dynamics. While our approach makes use of constrained optimization [31], the question 
at hand cannot be formulated as a simple optimization problem in terms of an aggregated 
objective function, such as the number of active nodes. Maximizing this number by ordi- 
nary means can lead to local minima or transient solutions that then fall back to asymptotic 
states with larger number of inactive nodes. A priori identification of the stable state that 
enjoys the desired properties is thus an important step in our formulation of the problem. 

Applications of our framework show that the approach is effective even when compen- 
satory perturbations are limited to a small subset of all nodes in the network, and when 
constraints forbid bringing the network directly to the target state. This is of outmost 
importance for applications because, in large networks, perturbed states in which the de- 
sired state is outside the range of admissible perturbations is expected to be the rule rather 
than the exception. From a network perspective, this leads to counterintuitive situations in 
which the compensatory perturbations oppose the direction of the target state and consist, 
for example, of suppressing nodes whose activities are already diminished relative to the 
target, but that lead the system to eventually evolve towards that target. These results 
are surprising in light of the usual interpretation that nodes represent "resources" of the 
network, to which we intentionally (albeit temporarily) inflict damage with a compensatory 
perturbation. From the state space perspective, the reason for the existence of such locally 
deleterious perturbations that have globally beneficial effects is that the target's basin of 
attraction, being nonlocal, can extend to the region of feasible perturbations even when the 
target itself does not. 

Our results based on randomly-generated networks also show that control is far more 
effective when focused on high-degree nodes. The effect is more pronounced in more het- 
erogenous networks, making it particularly relevant for applications to real systems, which 
are known to often exhibit fat-tailed degree distributions [32l [33] . More generally, we posit 
that the rate of successful control on a network will be strongly biased towards nodes with 
high centrality, including degree centrality, as shown here, and possibly other measures of 
centrality in the case of networks that deviate from random, such as closeness and between- 
ness. This suggests that in the situation where only a limited intervention is possible, one 
should focus primarily on central nodes in order to direct the entire network to a desired 
new state. We have demonstrated, moreover, that even when all nodes in the network are 
perturbed, the identification of compensatory perturbations is computationally inexpensive 
and scales well with the size of the network, allowing the analysis of very large networks. 

We have motivated our problem assuming that the network is away from its desirable 
equilibrium due to an external perturbation. In that context, our approach provides a 
methodology for real-time rescue of the network, bringing it to a desirable state before it 
reaches a state that can be temporarily or permanently irreversible, such as in the case 
of cascades of overload failures in power grids or cascades of extinctions in food webs, 
respectively. We suggest that this can be important for the conservation of ecological systems 
and for the creation of self-healing infrastructure systems. More generally, as illustrated in 
our genetic network examples, our approach also applies to change the state of the network 
from one stable state to another stable state, thus providing a mechanism for "network 
reprogramming" . 

As a context to interpret the significance of this application, consider the reprogramming 
of differentiated (somatic) cells from a given tissue to pluripotent stem cell state, which 
can then differentiate into cells of a different type of tissue. The seminal experiments 
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demonstrating this possibility involved continuous overexpression of specific genes [3^ . 
which is conceivable even under the hypothesis that cell differentiation is governed by the 
loss of stability of the stem cell state ^2]. However, the recent demonstration that the same 
can be achieved by temporary expression of few proteins [3S] indicates that the stem cell 
state may have remained stable (or metastable) after differentiation, allowing interpretation 
of this process in the context of the interventions considered here. This is relevant because 
a very small fraction of a population of treated cells is found to undergo the transition back 
to the stem cell state, suggesting that much could be learned from a modeling approach 
capable of identifying the perturbations most likely to direct the cellular network to adopt 
a new stable state. While induced pluripotency is an example par excellence of network 
reprograming, the same concept extends far beyond this particular system. Taken together, 
our results provide foundation for the control and rescue of network dynamics, and as 
such are expected to have implications for the development of smart traffic and power-grid 
networks, of ecosystems and Internet management strategies, and of new interventions to 
control the fate of living cells. 



Methods 

Constraints on incremental perturbations. The iterative procedure behaves well as 
long as the linearization involved in Eq. ^ remains valid at each step. The incremental 
perturbation at the point of closest approach, 5x(tc), selected under constraints (|2| alone 
will generally have a nonzero component along a stable subspace of the orbit x(t), which will 
result in 5xo larger than 5x(tc) by a factor 0(exp |A*'=(tc — i^o)|), where A*= is the finite-time 
Lyapunov exponent of the eigendirection corresponding to the smallest-magnitude eigen- 
value of M(xo,tc)- In a naive implementation of this algorithm, to keep ^xq small for the 
linear transformation to be valid, the size of 5x(tc) would be negligible, leading to negligible 
progress. This problem is avoided by optimizing the choice of 5x(tc) under the constraint 
that the size of (5xo is bounded above. Another potential problem is when the perturbation 
causes the orbit to cross an intermediate basin boundary before reaching the final basin of 
attraction. All such events can be detected by monitoring the difference between the linear 
approximation and the full numerical integration of the orbit, without requiring any prior in- 
formation about the basin boundaries. Boundary crossing is actually not a problem because 
the closest approach point is reset to the new side of the basin boundary. To assure that the 
method will continue to make progress, we solve the optimization problem under the con- 
straint that the size of 5xo is also bounded below, which means that we accept increments 
5xo that may temporarily increase the distance from the target (to avoid back-and-forth 
oscillations, we also require the inner product between two consecutive increments 5xo to 
be positive). These upper and lower bounds can be expressed as 

eo < |5xo| < ei. (7) 

This can lead to |(5x(tc)| ^ Ci due to components along the unstable subspace, but 
in such cases the vectors can be rescaled after the optimization. At each iteration, 
the problem of identifying a perturbation 5xo that incrementally moves the orbit toward 
the target under constraints and is then solved as a constrained optimization problem. 
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Nonlinear optimization. The optimization step of the iterative control procedure consists 
of finding the small perturbation (5xo that minimizes the remaining distance between the 
target, x*, and the system orbit x(t) at its time of closest approach, t^. Constraints are 
used to both define the admissible perturbations (|2| and also, as described in the previous 
paragraph, limit the magnitude of 5xo (j?]). The optimization problem to identify 5xo can 
then be succinctly written as: 



min 


|x*-(x(g + M(x'o,g -5x0)1 


(8) 


s.t. 


g(xo,Xo + 5xo) < 


(9) 




h(xo,Xo + 5xo) = 


(10) 




eo < |5xo| < ei 


(11) 




5xo • 5xq > 0, 


(12) 



where (12) is enforced starting from the second iteration, and 5xq denotes the incremental 



perturbation from the previous iteration. Formally, this is a nonlinear programming (NLP) 



problem, the solution of which is complicated by the nonconvexity of the constraint (11) 
(and possibly ^ and ([lo))). Nonetheless, a number of algorithms have been developed for 
the efficient solution of NLP problems, among them Sequential Quadratic Programming 
(SQP) [36j. This approach solves (|8| as the limit of a sequence of quadratic programming 
subproblems, in which the constraints are linearized in each sub-step. For all calculations, 
we used the SQP algorithm [37] implemented in the SciPy scientific programming package 
(http://www.scipy.org/). Note that this implementation does not require inverting matrix 
M (cf. Eq. 



Comparison with backward integration. Our approach should be compared with 
an apparently simpler alternative. Rather than keeping track of the variational matrix 
M(xo, t), which requires the integration of additional differential equations at every itera- 
tion, one could imagine using backwards integration of a trajectory starting at x(tc) + 5x(tc) 
to identify a suitable initial perturbation. This alternative procedure suffers the critical 
drawback that, for a particular choice of 5x(tc) (magnitude and direction), it is not certain 
that the time-reversed orbit will ever strike the feasible region defined by (|2]). This is 
particularly so in realistic situations where only a fraction of the nodes are accessible to per- 
turbations, resulting in a feasible region of measure zero in the full n-dimensional state space. 

Termination criteria and parameters. In all simulations, the control procedure is 
terminated if the updated initial condition attracts to within k = 0.01 of the target state 
within r = 10, 000 time units. Otherwise, we terminate the search if a compensatory 
perturbation is not found after a fixed number of iterations, which we set to be / = 1, 000. 
In general this number should be of the order of L/eQ, where L is the characteristic 
linear size of the feasible region. For each iteration, we use T = 10 time units within the 
integration step that identifies tc, which was estimated based on the time to approach the 
undesirable stable state. In the examples of Fig. [2} we used the parameters eo = 0.001 and 
ei = 0.01 for the optimization step. To lighten the computational burden in our Monte 
Carlo analysis, we used a more relaxed choice of eo = 0.005 and ei = 0.05 for the genetic 
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networks presented in Figs. |3]-[6| 

Construction of genetic networks. The networks are grown starting with a d-node 
connected seed network, by iteratively attaching a new node and connecting this node to 
each pre-existing node i with probabihty d x P^, where Yli^i ~ ^- '^^^ connections are 
assumed to be unweighted and undirected. We reject any iteration resulting in a degree-zero 
node, thereby ensuring that the final A^-node network is connected and has average degree 
> 2d for large N. Networks are generated for both uniform attachment probability, where 
Pi = l/N'"^ and A^** is the number of nodes in the network at iteration s, and for linear 
preferential attachment, where Pi = kf/Y2jkj and is the degree of node i at iteration 
s. The former leads to networks with asymptotically exponential distribution of degrees, 
which we refer to as homogeneous networks, while the latter leads to networks with asymp- 
totically power-law distribution of degrees, which we refer to as heterogeneous networks 
[32] . By construction, the resulting average degree is essentially the same for both homo- 
geneous and heterogeneous networks. In our simulations, we focused on networks with d = 2. 
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FIG. 1: Schematic illustration of the network control problem. (A) Network portrait. The goal is 
to drive the network to a desired state by perturbing a control set consisting of one or more nodes. 
(B) State space portrait. In the absence of control, the network at an initial state xq evolves 
to an undesirable equilibrium in the n-dimensional state space (red curve). By perturbing this 
state (orange arrow), the network reaches a new state Xq that evolves to the desired target state 
(blue curve). (C) Constraints. In general, there will be constraints on the types of compensatory 
perturbations that one can make. In this example, one can only perturb three out of n dimensions 
(equality constraints), which we assume to correspond to a thee-node control set, and the dynamical 
variable along each of these three dimensions can only be reduced (inequality constraints). This 
results in a set of eligible perturbations, which in this case forms a cube within the three-dimensional 
subspace of the control set. The network is controllable if and only if the corresponding slice of 
the target's basin of attraction (blue volume) intersects this region of eligible perturbations (grey 
volume). (D, E) Iterative construction of compensatory perturbations. (D) A perturbation to 
the current initial condition (magenta arrow) corresponds to a perturbation of the resulting orbit 
(green arrow) at the point of closest approach to the target. At every step, we seek to identify 
a perturbation that brings this point, and hence the orbit, closer to the target. (E) This process 
generates orbits that approach increasingly closer to the target (dashed curves), and is repeated 
until a perturbed state Xq is identified that evolves to the target. The resulting compensatory 
perturbation Xq — xq (orange arrow) brings the system to the attraction basin of the target without 
any a priori information about its location, and allows directing the network to a state that is not 
directly accessible by any eligible perturbation. 
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FIG. 2: Illustration of the control process in two dimensions. Yellow and blue represent the basins 
of attraction of the stable states xa and x^, respectively, while white corresponds to unbounded 
orbits. (A and B) Iterative construction of the perturbation for an initial state in the basin of 
with xs as a target (A), and for an initial state on the right side of both basins with x^ as 
a target (B). Dashed and continuous lines indicate the original and controlled orbits, respectively. 
Red arrows indicate the full compensatory perturbations. Individual iterations of the process are 
shown in the insets (for clarity, not all iterations are included). 
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FIG. 3: Control of genetic networks. (A) State space of the two-gene subnetwork represented 
in the inset, where the red curves mark the boundaries between the basins of x^, x^, and xc, 
and the arrows indicate the local vector field. (B) Illustration of compensatory perturbation on 
complex genetic networks, where each node is a copy of the two- gene system. We are given an 
initial network state xq representing the expression levels of the A'' gene pairs (color coded), and 
this state evolves to a stable state of the network x^ (top path). The goal is to knockdown one or 
more genes to reach a new state Xq that instead evolves to a target stable state x* 7^ x^ (bottom 
path). (C and D) Average computation time (C) and average number of iterations (D) required 
to control networks of N nodes with initial state xq = x^ and target x* = x^, demonstrating 
the good scalability of the algorithm. Each point represents an average over 1,000 independent 
realizations of homogeneous networks. 
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FIG. 4: Control of networks using a small fraction of all nodes. The dynamics follows ([6]) for e = 1 
and the control is illustrated for bringing the network from state to state x^. (A) Probability 
that an A^-node homogeneous network can be controlled by perturbing only k randomly selected 
nodes. (B) Corresponding results for heterogeneous networks with the same average degree (see 
Methods). Each point represents an average over 1,000 network realizations for one control set 
each. Thus, for control sets consisting of randomly selected nodes, heterogeneous networks are 
more likely to be controlled than homogeneous networks if the size of the control set is small. This 
should be contrasted with the case of control sets formed by the highest-degree nodes (see Fig. 
6), where heterogeneous networks are more likely to be controlled for both small and large control 
sets. 
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FIG. 5: Enrichment of high-degree nodes in successful control sets. (A) Probability that a given 
node of degree d is present in a successful control set for homogeneous (top) and heterogeneous 
(bottom) networks. This probability is scaled as P{d) = s{d)/M{d), where s is the number of times 
(counting multiplicity) a node of degree d appears in a successful control set, while M{d) is the 
total number of such nodes appearing in all successfully-controlled networks. (B) Dependence of 
the success rate on the average degree of the control set. Probability (red) for homogeneous (top) 
and heterogeneous (bottom) networks that a random set of 20% of the nodes can be used to control 
the corresponding network as a function of the average degree {d) in the set. The background of 
each panel (light blue) shows the distribution of {d) for randomly formed control sets. (C) Relation 
between success rate (left) and node degree (right) in an example network. The color bars show 
the percentage of successful control sets involving the given node (red) and the degree of the node 
(blue). Each point in panels A and B correspond to an average over 5,000 network realizations 
sampled 10 times while the statistics in panel C are based on sampling the given network 10, 000 
times, where = 50 and A: = 10 in all cases. The dynamics, coupling strength, and initial and 
target states are the same considered in Fig. |4j It follows that the ability to control a network 
correlates strongly with the degrees of the nodes in question. 
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FIG. 6: Increase in rate of successful control through high-degree nodes. (A and B) Probability of 
successfully controlling homogeneous (A) and heterogeneous (B) networks with a control set formed 
by a fraction of randomly selected nodes (gray) or by the same fraction of highest-degree nodes 
(orange). The fraction of nodes accessible in the control set is selected to be only 10% (circles) 
and 20% (diamonds) of the network. Each point represents an average over 2,000 independent 
network realizations for initial state and target state ic^, for the dynamics and coupling strength 
considered in Fig. |4j Biasing the intervention towards high-degree nodes dramatically improves 
the ability to control the network, permitting success with a small number of nodes where control 
through randomly selected nodes or low-degree nodes would fail. 
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